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HIGHLIGHTS 


►  The  OCFE  method  was  applied  to  efficiently  solve  the  diffusion  equation  in  a  solid  particle. 

►  The  OCFE  method  was  compared  to  some  other  methods  for  the  diffusion  equation  in  a  solid  particle. 

►  The  OCFE  method  was  used  to  solve  the  P2D  model  to  reduce  the  computation  time. 
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The  physics-based  pseudo  two-dimensional  (P2D)  model  for  lithium  ion  cells  requires  significant 
computation  time.  To  reduce  this  burden,  the  orthogonal  collocation  on  finite  elements  (OCFE)  method  is 
applied  here.  The  number  of  node  points  both  in  the  macro-scale  and  in  the  micro-scale  in  the  P2D 
model  is  reduced  by  optimally  locating  the  node  points.  This  OCFE  method  is  shown  to  be  better  than 
other  approximate  methods  for  the  spherical  diffusion  equation  with  flux  boundary  conditions.  This 
OCFE  method  is  also  used  to  solve  the  P2D  model  equations  and  the  results  are  compared  to  those 
obtained  from  COMSOL  (a  commercial  differential/algebraic  equation  solver  using  FEM). 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

A  model  with  accurate  short  time  response  is  important  for  the 
simulation  of  lithium  ion  cells  in  the  applications  with  high  rate 
charge/discharge  processes.  The  pseudo  two-dimensional  (P2D) 
model  developed  by  Doyle,  Fuller,  and  Newman  [  1,2  ]  has  been  widely 
used  to  simulate  the  performance  of  lithium  ion  cells  with  high  rates. 
The  P2D  model  is  a  physics-based  model  which  includes  material 
balances  and  charge  balances  in  the  solid  phases  and  the  electrolyte. 
The  governing  equations  for  the  concentration  of  the  electrolyte,  the 
solid  phase  potential,  and  the  electrolyte  potential  are  solved  in  the 
macro-scale  (the  direction  of  the  current  passing  through  the  cell). 
The  governing  equation  for  the  concentration  of  lithium  ions  in  the 
solid  phase  is  established  in  the  radial  direction  (the  pseudo  second 
dimension)  in  the  active  spherical  particles.  The  macro-scale  and  the 
micro-scale  (the  solid  particles)  are  coupled  through  the  equation  for 
the  transfer  of  lithium  ions  from  the  electrolyte  to  the  solid  particles. 
Both  the  macro-scale  and  the  micro-scale  equations  must  be  solved 
to  obtain  model  predictions.  Because  of  the  multi-scale  nature  of  the 
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P2D  model,  significant  computation  time  is  required  to  solve  the 
model  equations  using  numerical  methods  [3,4],  To  reduce  this 
computation  time,  several  approximation  methods  have  been 
applied  to  simplify  the  diffusion  equation  in  the  solid  phase.  These 
methods  include  the  polynomial  approximation  method  [5],  the 
residual  grouping  method  [6],  the  diffusion  length  method  [7],  the 
mixed  finite  difference  method  [8],  and  the  proper  orthogonal 
decomposition  (POD)  based  on  the  finite  volume  method  (FVM) 
[3,4],  In  this  paper  the  orthogonal  collocation  on  the  finite  elements 
(OCFE)  method  is  used  to  reduce  the  computation  time  by  reducing 
the  number  of  node  points  in  both  the  macro-scale  and  the  micro¬ 
scale  without  sacrificing  the  accuracy  of  the  predictions  from  the 
model.  These  methods  are  reviewed  for  spherical  diffusion  with  flux 
boundary  condition  and  compared  to  the  OCFE  method.  After  that 
the  OCFE  method  is  applied  to  the  complete  P2D  model  and  the 
results  are  compared  to  those  obtained  from  COMSOL 

2.  Transient  response  for  the  surface  concentration  of 
lithium  ions  in  a  spherical  particle 

The  governing  equation  of  the  concentration  of  lithium  ions  in 
a  spherical  particle  is  given  by  Fick’s  second  law: 
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where  Ds  is  the  diffusion  coefficient  of  lithium  ions  i 
At  the  center  of  the  particle,  there  is  no  molar  flux: 


Ds  ’  m  0 


The  molar  flux  at  the  surface  of  the  particle  equals  the  reaction 
rate: 
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cs(r,0)  =  0 
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The  difference  between  the  surface  concentration  and  the 
averaged  concentration  can  be  defined  as  follows: 

Acs,Rs(t)  =  cs(Rs,t)  -CSiavg(t)  (13) 

The  transfer  function  for  Eq.  (13)  is  given  by: 


AcSjrs  (s)  _  cs(Rs,s)  _  cs,avg(s) 

-JoT-^sT  i(s) 


(14) 


(3) 


where  Rs  is  the  radius  of  the  particle;  j(t)  denotes  the  current  density 
of  the  lithium  insertion/deintercalation  reaction;  as  is  the  specific 
area;  and  F  is  Faraday’s  constant.  Initially,  it  is  assumed  that  the 
concentration  of  lithium  ions  is  uniformly  distributed  in  the  particle: 

cs(r,0)  =  cSi  o 

First,  define  a  new  variable,  cs(r,  t)  =  cs(r,  t)  —  csq,  where  the 
initial  concentration  is  subtracted  from  the  concentration  in  the 
solid  particles.  The  model  equations  for  cs  are  similar  to  those  for  cs 
except  that  the  initial  condition  for  cs  is  zero,  and  are  presented  as 
follows: 


Note  that  Eq.  (14)  has  a  finite  steady  state.  Some  simplified 
models  for  spherical  diffusion  in  the  Laplace  domain  are  reviewed 
next. 

2.1.  Two-term  polynomial  approximation 

The  two-term  polynomial  approximation  [5]  to  Eqs.  (5)— (8)  is 
given  by  one  ordinary  differential  equation  and  one  algebraic 
equation: 


(4)  dCs,avg(t) 


Cs(Rs,  t)  -  Cs,avg(t)  =  -5DsQsP/'(f) 

The  transfer  function  for  Eq.  (16)  is: 

es(Rs,s)  -es,ayg(s)  Rs  n 

M  5DsasF  1 

2.2.  Three-term  polynomial  approximation 

The  three-term  polynomial  approximation  [5]  to  Eqs.  (5)— (8) 

d£s,avg(t) 


(15) 
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For  comparison  purpose,  the  transfer  function,  g(r,s),  of  cs(r,t) 
with  respect  to  the  input  j(t)  can  be  obtained  by  using  the  Laplace 
transform  method  as  follows  [9]: 


=  _ sinft(r^) _ 

asDsFr  Smh(Rsy/s/D?J  -Rsy/s/D~scosh(Rsy/s/D^ 

where  s  is  the  independent  variable  in  the  Laplace  domain.  The 
analytical  transfer  function  evaluated  on  the  particle  surface  is 
given  by: 


dq3VJt )  30DS  _  45 
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(18) 

(19) 


cs (Rs,  t )  -  cs,avg(t)  =  —RsqSVg(t)  -  p^pKt)  (20) 

where  qavg(t)  is  the  volume  averaged  flux  and  is  defined  in  Ref.  5 
with  the  initial  condition:  qaVg(0)  =  0.  The  transfer  function  for 
Eq.  (20)  is: 


cs(Rs,s)  —  Cs.avg(S) 
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Next,  define  the  averaged  value  of  cs(r,  t)  in  the  particle  as: 

Rs 

cs,avg(f)  =wf  r2^(r,t)dr  (11) 

s  o 

The  transfer  function  for  the  averaged  concentration  can  be 
obtained  by  substituting  Eq.  (9)  into  Eq.  (11),  that  is: 

Savg(s).  =  3  •  (12) 

j(s)  RsasFs  [1 


2.3.  Residues  grouping  [6]: 

The  analytical  transfer  function  for  the  surface  concentration 
difference  (Acs,Rs(t))  shown  in  Eq.  (14)  can  be  expanded  into  poles, 
Pic,  and  the  associated  residues,  resk,.  [6],  These  poles  and  residues 
are  sorted  such  that  the  values  of  the  poles  are  in  the  descending 
order.  The  analytical  transfer  function  shown  in  Eq.  (14)  can  be 
approximated  by  a  finite  number  of  poles  and  their  residues: 


^£sa(s)  ,  v^s- rt 

Us)  A^.s- 


resfc 


where  Z  is  the  steady  state  value  of  the  analytical  transfer  function; 
Pk  is  the  kth  pole  of  the  analytical  transfer  function;  and  resk  is  the 
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residue  associated  with  Pk.  The  expression  for  Z  can  be  obtained 
analytically: 


Z 


Rs 

"5  DsasF 


(23) 


The  N  poles  and  residues  in  Eq.  (22)  can  be  divided  into  groups 
(say  m  groups,  where  m  «  N).  In  this  work,  the  range  of  the  ten 
based  logarithm  for  the  absolute  values  of  the  N  poles  is  uniformly 
partitioned  into  m  segments.  The  corrected  residue  and  pole  for 
each  segment  are  given  by: 


reSgj 


resjf,  i  =  1-m 

(c  =  fcgi_1+l 


(24) 
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Hkikti  ,+l  PkTesk 
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where  kg,j,  i  —  0,1  is  the  index  for  the  bounds  of  the  pole 
segments.  Then,  the  approximation  to  Eq.  (22)  using  the  method  of 
residue  grouping  using  the  analytical  transfer  function  is  given  by: 


AcSiRs(s)  ,  ,  v'S-reSg* 
M  fc=is-Pg,k 


(26) 


2.4.  Finite  volume  method  [3]: 

The  geometry  discretization  for  a  particle  using  FVM  is  shown  in 
Fig.  1.  The  ID  geometry  is  divided  uniformly  into  n  control  volumes 
( CVs )  with  size  of  h.  The  node  points  are  located  at  the  center  of  the 
control  volumes.  Control  volume  “P”  is  represented  by  its  node 
point  “P”  and  has  two  faces:  face  “w”  and  face  “e”.  The  two  neighbor 
points  of  “P"  are  the  point  “VV”  and  the  point  “F\  It  is  assumed  that 
the  concentrations  in  the  control  volumes  are  piecewise  contin¬ 
uous,  that  is  the  concentration  in  control  volume  “P"  is  uniformly 
distributed  and  is  equal  to  the  concentration  at  node  point  “P”.  The 
discrete  model  using  the  FVM  is  obtained  by  integrating  Eq.  (5)  in 
each  control  volume  (see  Ref.  3  for  the  details).  The  resulting 
discrete  model  in  matrix  notation  is  given  by: 

C  =  AC(t)  +bj(t)  (27) 

where  C  =  [cs \ ,  cs 2,  Cs,n]T,  superscript  “T"  means  matrix 
transpose:  the  dot  over  C  means  time  derivative.  Matrix  A  in  Eq. 
(27)  is  an  n  x  n  tridiagonal  matrix  and  includes  the  diffusivity  and 
the  geometry  information.  The  concentration  on  the  particle 
surface  shown  by  an  open  circle  in  Fig.  1  at  r  =  Rs  is  not  solved 
directly  in  Eq.  (27)  but  it  can  be  obtained  by  applying  the  three- 


point  backward  difference  approximation  to  the  flux  on  the 
surface,  after  Eq.  (27)  has  been  solved: 

(28) 

The  Laplace  transform  of  Eq.  (27)  is  given  by: 


G(s)  =L.^'  =  (Is -A)-1  6  (29) 

where  G  (s)  is  a  n  by  1  vector  and  I  is  a  n  by  n  identity  matrix.  Based 
on  Eqs.  (12),  (28)  and  (29),  the  transfer  function  of  the  surface 
concentration  difference  using  the  FVM  is  given  by: 


Acs,Rs(s) 

:I(S) 


3  h 

8  DsasF 


3 

RsasFs 


(30) 


where  Gn-t(s)  and  G„(s)  are  the  (n-l)th  and  the  nth  element  of  G(s) 
inEq.  (29). 


2.5.  Residue  grouping  using  the  state  space  approach  [6]: 


The  transfer  function  shown  in  Eq.  (29)  has  a  zero  pole  (s  =  0).  To 
get  rid  of  this  pole  at  the  origin,  the  first  equation  was  subtracted 
from  Eq.  (2)  through  n  shown  in  Eq.  (27).  The  resulting  equations 
are  given  by: 

X  =  A  x(t)  +  bj(t)  (31) 

where  =  cs  l+1  -  cs  l .  i  —  1, ...,  n  -  1  ;  A  is  an  n-1  by  n-1  matrix 
with  A jj  =  Ai+lj+1  -  A]  J+1 ,  i.j  —  1. ....  n  -  1  ;  and  b  is  an  n-1  by  1 
vector  with  bj  =  bi+1  —  b1 ,  i  =  l,...,n  — 1.  Next  define  the 
output  as: 

l(t)  =  C  x(t)  (32) 

where  C  is  an  n-1  by  n-1  identity  matrix.  The  transfer  function  for 
this  state  space  approach  is  given  by: 


y(s) 

j(S ) 


where: 


l  = 


-CA  'b 


(33) 


(34) 


rk=- k=  1,2,..., n-1  (35) 

qk  is  the  /cth  eigenvector  associated  with  the  kth  eigenvalue  Xu  in 
the  expansion  of  A  : 


Fig.  1.  Discretized  geometry  for  the  diffusion  problem  using  the  finite  volume  method 
(It  =  0.5  h,  12  =  1.5h). 


AQ.  =  Q4  (36) 

q k  is  the  kth  column  in  Q;  A  is  a  diagonal  matrix  and  its  non-zero 
elements  are  the  eigenvalues  of  matrix  A.  The  vector  pk  is  the  kth 
row  in  matrix  P  where  P  =  Q-1.  In  the  modal  form,  the  matrices  A, 
b,  and  C  in  Eqs.  (31)  and  (32)  are  equivalent  to  : 


A*  =  diag(A],  A2,...,An_1)  =  A 

(37) 

b  =  [1,  1  -,l]T 

(38) 

C  =  {lift,  X2r2,...,Xn_-lr„_i] 

(39) 
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Next,  sort  the  eigenvalues  of  A  in  decreasing  order  and  change 
the  sequence  of  the  vectors  in  Q  and  pk  in  P  correspondingly. 
Evaluate  the  residues  according  to  Eq.  (35)  and  update  Eq. 
(37)-(39).  Divide  the  eigenvalues  into  Nrgss  segments.  In  a  manner 
similar  to  the  method  of  the  residue  grouping  with  the  analytical 
transfer  function,  the  grouped  residues  and  eigenvalues  are  given 
as  follows: 

ksi 

Lg,i  =  £  Lk,  i=  1,...  ,NRGss  (40) 

k = kEjj  i+1 

Ag.  Vk,n-1^  f _  r  Nrgss  (41) 

The  approximate  transfer  function  using  the  residue  grouping 
with  state  space  (RGSS)  approach  is  given  by: 


l42) 

where  Zg  is  the  approximation  of  Z  obtained  by  substituting  Eqs. 
(40),  (41 )  to  Eq.  (37)— (39).  Recall  that  the  dependent  variable  in  Eq. 
(31 )  is  the  difference  between  the  deviation  concentration  at  the 
node  points  except  for  the  first  point  (csj,  i  =  2 ,...,n)  and  the 
deviation  concentration  at  the  first  node  point  (cs  l ).  The  transfer 
function  of  the  deviation  concentration  at  the  first  node  point 
obtained  from  the  previous  section  (section  of  Finite  Volume 
Method)  should  be  added  to  the  transfer  function  shown  in  Eq.  (42) 
to  obtain  the  transfer  function  for  each  of  the  concentration  devi¬ 
ations  at  the  node  points  from  2  to  n.  Similarly,  the  transfer  function 
of  the  surface  concentration  difference  can  be  obtained  by  using 
a  three-point  backward  difference: 

^s,Rs(s)  9  1  3  h  ,  3 

1(ip  ~  8Grgss  "  1  (S)  “  8Grgss’"-2(S)  “  8Drf  +  W 


The  frequency  response  of  the  surface  concentration  difference 
obtained  by  using  the  various  models  shown  above  can  be  calcu¬ 
lated  by  replacing  Laplace  s  by  ito  where  i  is  the  imaginary  unit  and 
<,)  denotes  the  frequency  in  rad  s-1  Fig.  2  shows  a  comparison  of  the 
magnitudes  of  the  frequency  response  between  the  analytical 
solution  and  the  simplified  models.  There  are  an  infinite  number 
of  negative  poles,  which  are  the  roots  of  the  eigen-equation:  tan 
/i(Rss/Ds)  =  Rss/Ds  in  the  analytical  transfer  function  of  the  surface 
concentration  deviation  shown  in  Eq.  (10).  The  solid  line  shown  in 
Fig.  2  denotes  the  analytical  solution  of  the  frequency  response.  The 
dashed  line  shown  in  Fig.  2  denotes  the  results  of  the  truncated 
analytical  solution  with  the  first  50  slowest  (smallest  in  magnitude) 
poles.  For  the  long  time  scale  (small  values  of  w)  the  magnitude  of 
the  analytical  frequency  response  approaches  to  a  stable  value  as 
expected.  As  the  frequency  increases,  the  logarithm  of  magnitude 
of  the  analytical  frequency  response  goes  to  a  very  small  value.  The 
boxes  represent  the  frequency  response  of  the  9  terms  grouping  of 
the  50th  order  truncation  of  the  analytical  solution.  The  frequency 
response  obtained  using  this  residue  grouping  method  with  9 
groups  for  the  50th  order  truncation  agrees  well  with  the  analytical 
solution  when  the  frequency  is  less  than  about  10  rad  sec-1.  Due  to 
the  loss  of  the  fast  poles,  the  frequency  response  of  the  50th  order 
truncation  of  the  analytical  solution  tends  to  a  constant  value  when 
the  frequency  is  greater  than  10  rad  sec-1.  The  frequency  responses 
of  the  finite  volume  method  with  30,  50,  and  80  control  volumes 
are  represented  by  stars,  a  dotted  line,  and  crosses,  respectively  in 


Frequency  a>,  rad  sec'1 


Frequency  response  of  the  simplified  models  compared 


the  analytical 


Fig.  2.  The  frequency  responses  of  all  of  these  FVM  cases  overlap  the 
analytical  solution  when  the  frequency  is  less  than  1  rad  sec-1.  As 
the  number  of  control  volumes  increases,  the  accuracy  of  the 
frequency  response  of  the  FVM  in  the  short  time  scale  (high 
frequency)  improves.  But  the  FVM  with  80  control  volumes  is  less 
accurate  than  the  50th  order  truncation  of  the  analytical  solution  in 
the  short  time  response  («  >  10  rad  sec-1).  The  circles  shown  in 
Fig.  2  denote  the  frequency  response  of  the  residue  grouping  using 
the  state  space  approach  (RGSS).  The  state  space  was  obtained  by 
using  FVM  with  50  CVs.  The  frequency  response  of  the  RGSS  is  close 
to  that  of  the  FVM  with  50  CVs  when  the  frequency  is  less  than 
0.01  rad  sec-1  or  the  frequency  is  greater  than  3  rad  sec-1.  In  the 
frequency  range  [0.01  3]  rad  sec-1,  the  frequency  response  of  FVM 
agrees  well  with  the  analytical  solution,  but  the  frequency  response 
of  the  RGSS  oscillates  around  the  analytical  solution.  The  frequency 
responses  of  the  two-term  polynomial  and  the  three-term  poly¬ 
nomial  approximations  are  denoted  by  the  up  triangles  and  dia¬ 
monds  in  Fig.  2,  respectively.  As  shown  in  Fig.  2,  the  two-term 
polynomial  approximation  only  captures  the  long  time  steady 
state  response.  The  three-term  polynomial  approximation  is 
improved  somewhat  by  adding  one  pole  to  the  two-term  poly¬ 
nomial  approximation.  The  pole  is  given  by  -30 Ds/Rs  [2]  which 
equals  to  6.0  x  10-3  s-1  based  on  the  parameters  shown  in  Table  1. 
This  pole  is  between  the  first  and  the  second  analytical  poles  which 
were  obtained  by  solving  the  eigen-equation. 

3.  Orthogonal  collocation  on  finite  elements 

To  use  the  orthogonal  collocation  on  finite  elements  technique, 
the  geometry  (radial  direction  in  the  solid  particle)  is  divided  into 


Table  1 

Parameters  for  solid  phase  diffusion  [6], 
Parameter 

Diffusion  Coefficient  Ds,  m2  s-1 
Radius  of  the  Particle  Rs,  m 
Specific  Area  as,  m2  m  3 
Faraday's  Constant  F,  C  mor1 
C  Rate  Current,  A 


2.0  x  10~16 
1.0  x  10-6 
1.74  x  106 
96,487 
0.272 
1.18 
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number  of  elements  shown  in  Fig.  3.  In  each  element,  the  colloca¬ 
tion  points  (i.e.,  the  inner  points  in  an  element)  are  located  at  the 
roots  of  the  Jacobi  polynomial  whose  order  equals  to  the  number  of 
collocation  points  in  the  element.  The  IVth  order  Jacobi  polynomial 
is  defined  by  [10]: 


(,  -  xjVpT,  .  [<*  -  xTVWJ ,  X.  [0, 1] 

(44) 

where  a  and  0,  coefficients  of  Jacobi  polynomial,  determine  the 
locations  of  the  collocation  points.  When  the  value  of  a  is  increased, 
the  locations  of  the  inner  points  in  each  element  move  closer  to  the 
left  boundary  of  the  element  (eg  in  Fig.  3).  When  the  value  of  /?  is 
increased,  the  locations  of  the  inner  points  move  closer  to  the  right 
boundary  of  the  element  (rf  in  Fig.  3).  When  both  a  and  /3  equal 
zero,  the  inner  node  points  are  symmetric  about  the  center  of  the 
element  as  shown  in  Fig.  3.  Each  element  is  rescaled  such  that  the 
local  coordinate  in  the  element  ranges  from  0  to  1.  As  shown  in 
Fig.  3,  the  element  “e”  with  the  size  of  he  starts  at  rg  and  ends  at  rf. 
The  relationship  between  the  local  coordinate  ?  and  the  global 
coordinate  r  is  given  by: 


d2cf 
dr 2 


(49) 


The  locations  of  the  node  points  can  be  determined  by: 


r(e  l)(Noc:1)tj  —  r0  +  ke^OCj  (50) 

where  e  =  1 . NE,  j  =  2,..., Noc+2,  n  =  0,  R0c  is  a  vector  of  the 

locations  of  the  roots  of  Noe  order  Jacobi  polynomial  including  the 
two  boundaries  of  the  element.  The  Noc+2  by  Noc+2  matrices  Aoc, 
Boc  and  the  Noc+2  by  1  vector  Roc  can  be  determined  by  invoking 
the  routines  in  Ref.  10  (i.e.,  JCOB1  presented  in  Appendix  A1  in 
Ref. 10).  The  model  equations  at  the  inner  points  (collocation  points) 
in  each  element  were  obtained  by  substituting  Eq.  (48)-(50)  to  Eqs. 
(5)— (7).  In  the  interfaces  between  the  two  elements,  the  mass  flux 
is  continuous.  The  summary  of  the  model  equations  using  the 
orthogonal  collocation  on  finite  elements  is  given  as  follows. 

At  the  ith  collocation  point  for  element  e: 


deg. 

dt 


DsB0c(t,:)-cf  +  r - 2Ds - AoeftO-cg  i 

r(e— 1)(JV0C+1)+! 

2,...,N0C  + 1 


(51) 


The  concentration  at  the  node  points  in  the  element  e  is 
denoted  by  cfj  where  j  =  1,..., Noc+2,  and  Noc  is  the  number  of 
inner  points.  If  the  same  number  of  inner  node  points  are  assigned 
to  all  the  elements,  the  total  number  of  node  points  in  the  whole 
geometry  is:  NT  =  NE  x  (N0c+1)+1,  where  NE  denotes  the  total 
number  of  elements.  Globally,  the  concentration  at  the  node  points 
is  denoted  by  Q  where  i  =  1,..,Nt.  For  element  e  with  Noc  =  4, 
the  first  derivative  and  the  second  derivative  of  the  local  concen¬ 
tration  with  respect  to  the  local  coordinate  can  be  approximated 
by:  [10]: 


where  AocO'.O  is  a  row  vector,  and  i  denotes  the  ith  row  in  matrix 
Aoc,  cf  is  the  Noc+2  by  1  column  vector  for  the  local  concentration 
and  the  dot  between  two  vectors  denotes  the  dot  product. 

At  r  =  rf  for  any  element  beside  of  the  element  NE: 

-^AociNoc  +  2,:)  -cf  =  -^AocJl , :)  -cf 1  (52) 

he  he+1 

At  r  =  0  (r  =  r0e  and  e  =  1): 

-J£aoc(1,:).cJ  =  0  (53) 

At  r  =  Rs  (r  =  rle  and  e  =  NE): 


dg 

dt; 


=  Aoccf 


(46) 


dt 


(47) 


where  cf  =  [cf  1,cf2,cf3,cf4,cf5,cf6]T,  cf1  andcf6  represent  the 


concentrations  at  rg  and  rf ,  respectively.  Converting  the  derivatives 
for  the  concentration  from  the  local  coordinate  to  the  global  coor¬ 
dinate  according  to  the  relationship  between  the  local  coordinate  to 
the  global  coordinate  shown  in  Eq.  (45),  we  have: 


d<l 

Hr 


(48) 


Cj  Ci+6 

,  c',f  he  ,  , 

— *  »j»  * — »  •  t •  • — * 

0  r£ _ ?  rl  Rs 


Fig.  3.  Schematic  of  the  node  points  assignment  for  the  diffusion  problem  using  the 
method  of  orthogonal  collocation  on  finite  elements  (number  of  elements:  NE  =  4; 
number  of  inner  points:  N0c  =  4;  a  =  0  and  0  =  0). 


-^Aoc(Noc  +  2,:)  .&=fF  (54) 

The  discrete  model  using  OCFE  shown  in  Eq.  (51)-(54)  for  the 
1-D  diffusion  problem  can  be  rewritten  as: 


MC  =  AC+bj(t)  (55) 

where  the  mass  matrix  M  is  a  modified  identity  matrix  in  which  the 
diagonal  element  is  1  for  the  collocation  points  and  is  zero  on  the 
interfaces  between  the  elements.  The  coefficient  matrix  A  in  Eq. 
(55)  is  a  banded  matrix  with  the  bandwidth  of  2(Noc+l  )+l.  The 
Laplace  transform  of  Eq.  (55)  is  given  by: 

CocfeW  -  II  =  (sM  +  ArVb  (56) 

The  transfer  function  of  the  difference  between  the  surface 
concentration  and  the  average  concentration  using  the  OCFE 
method  is  given  by: 


Acs,Rs(s) 

m 


=  Gqcfe,n,(s)  + 


3 

RsasFs 


(57) 


The  frequency  response  of  Acs  Rs  using  OCFE  can  be  obtained  by 
substituting  iw  for  s.  Fig.  4  shows  the  frequency  response  of  AcSiRs 
using  OCFE  with  uniform  elements,  a  different  number  of  elements, 
and  a  different  number  of  inner  points.  For  the  first  set  of  cases, 
the  coefficients  a  and  fl  used  in  the  OCFE  were  set  to  zero.  The 
solutions  of  OCFE  were  compared  to  the  solution  of  FVM  with  30 
control  volumes  and  the  analytical  solution.  The  numbers  in  the 
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Frequency  £u,  rad  sec'1 


Fig.  4.  Effects  of  the  number  of  elements,  the  number  of  collocation  points,  the  value 
of  jS,  and  the  size  of  the  element  near  the  particle  surface  on  the  frequency  response  of 
the  surface  concentration  using  OCFE. 


parenthesis  shown  after  the  model  names  in  Fig.  4  indicate  the 
number  of  equations  in  the  corresponding  discrete  models.  With 
a  similar  number  of  equations,  the  discrete  model  using  OCFE  with 
more  inner  points  possesses  higher  accuracy  in  the  frequency 
response  of  Acs,rs.  When  the  number  of  inner  points  is  equal  to  4, 
the  OCFE  model  with  21  equations  (circles  in  Fig.  4)  is  more  accu¬ 
rate  than  the  FVM  with  30  equations  (filled  triangles  in  Fig.  4)  in  the 
short  time  response  (oj  =  1  rad  sec-1)  of  the  surface  concentration 
relative  to  the  average  concentration. 

Fig.  4  also  shows  the  effect  of  the  value  of  (3  on  the  accuracy  of 
the  short  time  response  of  the  surface  concentration.  Four  equal 
sized  elements  were  assigned  to  the  geometry.  In  each  element, 
four  collocation  points  were  applied  and  the  value  of  a  was  fixed 
to  zero.  The  frequency  responses  using  the  OCFE  with  four 
different  values  for  /?:  0,  2,  5,  and  10  were  compared  to  the 
analytical  solution.  As  the  value  of  ft  increases,  the  short  time 
response  of  the  surface  concentration  relative  to  the  average 
concentration  gets  closer  to  the  analytical  solution.  This  is 
because  the  higher  the  value  of  fl  is,  the  closer  the  collocation 
points  are  to  the  right  boundary  of  the  element.  For  the  element 
near  the  particle  surface,  high  value  of  /?  leads  to  the  locations  of 
the  collocation  points  moving  to  the  surface  of  the  particle,  which 
improves  the  accuracy  of  the  approximation  of  the  gradient  at 
the  surface  shown  in  Eq.  (54).  Smith  and  Wang  [11  indicated  that 
placing  a  smaller  element  near  the  particle  surface  helps  to 
achieve  a  better  approximation  to  the  surface  concentration.  In 
this  work,  we  used  unequal  sized  elements  for  the  geometry  to 
obtain  more  accurate  short  time  (large  w  values)  response  of  the 
surface  concentration  relative  to  the  average  concentration.  To  do 
this,  we  defined  a  factor,  /  such  that  the  size  of  the  element  near 
the  surface  is  determined  by: 

hNE  =  (1  -f)Rs  (58) 

The  size  of  each  of  the  other  elements  is  identical  and  is  deter¬ 
mined  by: 

he  =jv^TT  e  =  V->ne-1  (59) 

Fig.  4  shows  the  effects  of  the  factor/on  the  frequency  response 
of  the  surface  concentration  relative  to  the  average  concentration. 


The  values  of  the  coefficients  a  and  (3  were  set  to  zero.  The  geometry 
was  partitioned  into  four  elements,  each  of  which  includes  four 
collocation  points.  The  frequency  responses  using  the  OCFE  with 
different  values  of  /  were  compared  to  the  analytical  solution  in 
Fig.  4.  The  factor  /  is  equal  to  0.75  for  the  case  using  uniform 
elements  with  Ng  —  4.  When /is  increased  the  size  of  the  element 
near  the  surface  decreases,  and  the  accuracy  of  the  short  time 
response  of  the  surface  concentration  relative  to  the  average 
concentration  increases.  When/is  equal  to  or  greater  than  0.96,  the 
accuracy  of  the  short  time  response  of  the  surface  concentration 
relative  to  the  average  concentration  is  improved  significantly. 
Fig.  4  also  shows  that  when  /  =  0.99  the  response  of  the  surface 
concentration  relative  to  the  average  concentration  is  much  more 
accurate  than  that  when  /  =  0.96  for  w>=5  rad  sec-1.  But  the 
response  of  the  surface  concentration  relative  to  the  average 
concentration  is  worse  in  the  frequency  region  around  1  rad  sec-1 
when  /  =  0.99  compared  to  the  case  with  f  =  0.96.  This  can  be 
corrected  by  increasing  the  number  of  elements  or  increasing  the 
number  of  inner  points  in  each  element. 

The  method  of  OCFE  was  applied  to  solve  the  P2D  model  for 
a  lithium  ion  cell  with  mesocarbon  microbead  (MCMB)  as  the 
anode  active  material  and  lithium  cobalt  oxide  as  the  cathode 
active  material,  and  1M  LiPF6  in  propylene  carbonate/ethylene 
carbonate/dimethyl  carbonate  mixture  as  the  electrolyte.  The 
parameters  for  the  P2D  model  which  were  presented  in  our 
previous  work  [3]  and  are  summarized  in  Table  2.  The  open  circuit 
potentials  of  the  active  material  in  the  electrodes  were  given  in 
Ref. 14.  The  expressions  of  the  diffusion  coefficient  of  lithium  ions  in 
the  electrolyte  and  the  ionic  conductivity  of  the  electrolyte  were 
presented  in  Ref. 15.  The  particles  in  the  electrodes  were  divided 
into  four  elements  in  each  of  which  four  collocation  points  were 
assigned  with  a  =  0  and  /?  =  0.  The  factor /was  set  to  0.96  for  the 
discretization  of  the  radius  of  the  particle.  The  positive  electrode 
and  the  negative  electrode  were  partitioned  into  five  equal  sized 
elements.  And  one  element  was  assigned  to  the  separator.  In  each 
element  in  the  both  electrodes  and  the  separator,  five  collocation 
points  were  assigned  with  a  =  0  and  /}  =  0.  The  resulting  discrete 


Table  2 

Parameters  for  the  P2D  model. 


Parameter  Anode  (LixCf0 


Thickness,  pm 
Particle  radius,  pm 
Solid  phase  diffusion 
coefficient,  m2  s  ' 

Reaction  rate  constant, 

Maximum  concentration  in 
solid  phase,  mol  m  3 
Porosity 

Volume  fraction  of  fillers 
Electronic  conductivity,  S  m-1 
Charge  transfer  coefficient 
Bruggeman’s  number 
Initial  state  of  charge 
Temperature,  ”C 
Thickness  of  separator,  qm 
Porosity  of  separator 
Bruggeman's  number  in  separator 
Cationic  transference  number 

electrolyte,  mol  nr3 
C  rate  current  density,  A  cm~2 
a  Obtained  from  Ref.  12. 
b  Obtained  from  Ref.  13. 
c  Obtained  from  Ref.  14. 


73.5 

12.5 

5.5  x  io-,4a 
1.764  x  10-llb 
30556c 

0.4382b 

0.0566b 

0.5a 

4.1b 

0.756 

25 

25b 

0.45b 

2.3b 

0.435b 

1000 

25 


Cathode 

(LiyCoOj) 

"700 

8.5 

1.0  x  10-1Ia 
6.667  x  10-nb 


51555c 

0.3b 
0.1 5b 
10a 
0.5a 
1.5b 
0.465 
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rates  varying  from  C/10  to  20C  (lines  denote  the  results  from  COMSOL  Mul- 
;  and  symbols  are  the  results  obtained  by  using  the  OCFE  method). 


model  using  OCFE  involving  1608  equations  were  solved  using  the 
DAE  solver  DASRT  under  different  C  rates  ranging  from  C/10  to  20C. 
The  P2D  model  was  also  solved  in  COMSOL  Multiphysics®  [16], 
a  commercial  finite  element  solver,  with  the  degrees  of  freedom  of 
11,270  (number  of  unknown  variables  to  be  solved)  and  under  the 
same  conditions  as  in  the  OCFE  model.  The  discharge  profiles  for 
both  models  under  the  selected  C  rates  are  shown  in  Fig.  5.  The 
symbols  denote  the  results  of  the  OCFE  model  and  the  lines  denote 
the  results  of  COMSOL  Multiphysics®  model.  The  two  models  agree 
well  in  the  current  rate  range:  C/10-20C.  And  the  OCFE  model  is 
about  30  times  faster  than  the  COSMOL  Multiphysics  model 
(approximately  2  s  vs.  60  s  computation  time). 

4.  Conclusions 

In  this  paper,  the  methods,  which  include  the  polynomial 
approximation,  the  residue  grouping,  and  the  finite  volume 
method,  to  solve  the  1-D  diffusion  problem  for  the  frequency 
response  of  the  surface  concentration  of  the  particle  were 
compared  to  the  analytical  solution  for  spherical  diffusion  with  flux 
boundary  conditions.  According  to  the  comparisons  of  the 
frequency  response  of  the  surface  concentration  relative  to  the 
average  concentration,  it  was  found  that  the  two-term  polynomial 
approximation  only  captures  the  steady  state  of  system  and  the 
three-term  polynomial  approximation  provided  more  accuracy  by 
adding  one  pole  to  the  transfer  function  of  the  surface  concentra¬ 
tion.  The  residue  grouping  with  the  analytical  transfer  function 
approach  agrees  well  with  the  high  order  truncation  of  the 
analytical  transfer  function.  The  accuracy  of  the  residue  grouping 
with  the  state  space  approach  depends  on  the  accuracy  of  the 
original  finite  volume  method  upon  which  the  state  space  model 
was  developed. 

The  OCFE  method  was  introduced  in  this  paper  to  solve  the 
diffusion  problem  in  the  solid  particle.  It  was  found  that  when  the 
total  number  of  node  points  is  fixed,  increasing  the  number  of 
collocation  points  by  decreasing  the  number  of  elements  improves 
the  accuracy  of  the  approximation  of  the  surface  concentration  in 
the  short  time  scale.  It  was  also  found  that  reducing  the  size  of  the 
element  near  the  surface  of  the  particle  significantly  improves 
the  accuracy  of  the  approximation  of  the  surface  concentration  in 


the  short  time  period.  The  OCFE  method  was  also  successfully 
applied  to  solve  the  P2D  model  for  a  lithium  ion  cell  and  reduced 
the  computation  time  by  approximately  30  times  compared  to 
COMSOL  Multiphysics®. 
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List  of  symbols 

as  pecific  area  of  the  electrode,  m2  m  3 

A  coefficient  matrix 

Aoc  matrix  in  the  approximation  of  the  first  order  derivative 
using  OCFE  method 
b  source  term  vector 

Boc  matrix  in  the  approximation  of  the  second  order 

derivative  using  OCFE  method 

cs  concentration  of  lithium  ions  in  solid  phase,  mol  m  3 

C  vector  of  concentration  at  the  node  points,  mol  m-3 

Ds  diffusion  coefficient  of  lithium  ions  in  solid  phase,  m2  s  1 

F  Faraday’s  constant,  96,487  C  mol-1 

G  transfer  function  vector 

h  size  of  the  control  volume 

he  size  of  the  element  e 

1  identity  matrix 

fapp  applied  current,  A 

j(t)  time  dependent  volumetric  current  density,  A  m-3 

M  mass  matrix 

Ne  number  of  elements 

Noc  number  of  collocation  points 

Nj  total  number  of  node  points 

p  pole  of  the  transfer  function 

Roc  vector  of  the  positions  of  the  collocation  points 

Rs  radius  of  the  particle,  m 

s  independent  in  Laplace  domain 

Z  steady  state  of  the  transfer  function 

Greek 

a  coefficient  of  Jacobi  polynomial 

/3  coefficient  of  Jacobi  polynomial 

?  local  coordinate  in  radial  direction  in  the  particle 

w  frequency,  rad  s-1 

Subscript 

avg  average 

Rs  evaluated  at  the  surface  of  the  particle 
max  maximum 

Superscript 

-1  matrix  inverse 

e  element  index 

T  matrix  transpose 
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